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Much recent discussion about dark matter has been centered on two seemingly independent prob- 
lems: the abundance of substructure in dark matter halos, and the cuspiness of the halos' inner 
density profile. We explore possible connections between the two problems by studying the gravi- 
tational scattering effects due to subhalos on the phase-space distribution of dark matter particles 
in the main halos. Our series of controlled numerical experiments indicates that the number and 
mass density of subhalos can be high enough to cause the collisionless dark matter particles in the 
inner part of a main halo to diffuse, flattening the main halo's inner cusp within a few dynamical 
times. Depending on the masses and concentration of the subhalos, the inner density profile of the 
whole system (main plus sub halos) can either steepen or flatten. Subhalo accretion can therefore 
introduce significant scatter in the inner density profiles of dark matter halos, offering a possible 
explanation for the range of profiles seen in both observations and cosmological simulations. 



Particles undergo random walks and diffusion through 
collisional scatterings. The most noted example is the 
Brownian motion of small macroscopic particles, whose 
velocities exhibit frequent sudden changes due to impul- 
sive collisions with individual molecules in a liquid. On 
astrophysical scales, stars also undergo random walks in 
velocity space due to gravitational scatterings with, e.g., 
other stars or giant interstellar clouds in a galaxy |lj. 

Recent high resolution TV-body simulations of hierar- 
chical structure formation in cold dark matter (CDM) 
models have shown that spatial distribution of dark mat- 
ter in galaxy-hosting halos is not entirely smooth. In- 
stead, roughly 10% of a halo's mass is in the form of hun- 
dreds to thousands of smaller, dense satellite subhalos of 
varying mass [2j ■ In this Letter we examine whether these 
subhalos can be the source of a fluctuating gravitational 
potential that produces collisional transport of CDM par- 
ticles in the main halo, even when the self-interaction of 
CDM is collisionless. Our approach is based on numeri- 
cal simulations and addresses the fully nonlinear regime 
of halo-subhalo interaction; a complementary approach is 
pursued by |3| , who have used second-order cosmological 
perturbation theory to derive a kinetic equation for the 
phase-space distribution of halo dark matter particles. 

Physics of Diffusion. — A test particle of mass M t and 
velocity Vt experiences dynamical friction and exhibits 
random walks (in velocity space) as it moves through 
the gravitational potential of background particles of 
mass Mft. Both processes change the test particle ve- 
locity (Avi,i — 1,2,3) and energy (AE): the dynamical 
friction is described by the diffusion coefficient D(Avii) 
where dv t /dt = v t D(Av\\); the random walk is described 
by the diffusion tensor D(AviAvj). The rate of change 
of the kinetic energy of the test particle is D(AE) = 
M t Hi[vi D(Avi) + \D(Av 2 )] 0. For background par- 
ticles with a uniform mass density pb and an isotropic 
Maxwellian velocity distribution with dispersion o~b, it is 

D(AE) = 4TrG 2 ^^\nA{-M t F{x) + M b [er£(x) - F(x)]} 

(1) 



where In A is the Coulomb logarithm, x = v t /V%o~b, and 
F(x) — erf (a;) — 2x exp(— x 2 )/y/ir. Note that this equa- 
tion is valid for arbitrary ratios of Mt/Mb 0- The first 
term in Eq. describes the energy loss of the test par- 
ticle due to dynamical friction. In the standard Chan- 
drasekhar picture, a large test mass M t scatters off a sea 
of small background particles with mass Mf,. In this limit 
(M t S> Mb), the first term in Eq. (due to the test par- 
ticle polarizing the background medium) dominates, and 
the second term is typically ignored. 

Our focus in this Letter is different. We are interested 
in the effects on the dark matter particles in the main 
halo (our test particles) due to the ensemble of subha- 
los (our background particles). The relevant mass range, 
Mt <C Mb, is therefore opposite of that in the last para- 
graph. The polarization cloud term is completely neg- 
ligible. Instead, the key process is the second term in 
Eq. (JTJ , which describes the heating of the test particle 
due to stochastic fluctuations in the background parti- 
cles. Changes in the potential due to the distribution 
of dark matter substructure are the dominant scattering 
source in our study. 

Effects of Diffusion. — We perform a series of fully dy- 
namical simulations using GADGET, a publicly available N- 
body tree code , to follow the evolution of dark matter 
in a parent halo containing an ensemble of subhalos. To 
study the dynamical interplay between a main halo and 
its subhalos in a controlled and semi-realistic way, we use 
subhalo properties similar to those from earlier full-scale 
cosmological simulations 0. This strategy allows us to 
perform a suite of numerical experiments to quantify the 
effects due to a wider range of subhalo masses, concentra- 
tion, and orbits than is possible with large cosmological 
simulations. 

Initially the particles in the main halo are given an 
NFW radial density profile 0|: p(r) = p crtt S/[x(l + x) 2 ], 
where x = r/r s ,S = 200c 3 /3[ln(l + c) - c/(l + c)], and 
the concentration parameter c = r V i r /r s is the ratio of 
'the halo's virial to scale radius. We use a total of 10 6 
simulation particles for the main halo and a force soften- 



2 



ing of 0.0157V The particle velocities are drawn from a 
local isotropic Maxwellian distribution where the radius- 
dependent velocity dispersion is computed from the Jeans 
equation (see @ for details and tests of a similar set-up). 
This velocity setup may cause an NFW halo's inner cusp 
to artificially flatten initially 9J. To work around this 
problem, we evolve an NFW halo in isolation for ~ 8 td yn 
(where t 2 dyn = 3n/16Gp(r s )) and use this evolved halo as 
our initial main halo. Our tests have shown that this 
halo, which does differ slightly at r < 0.lr s from its orig- 
inal NFW structure, is extremely stable over the next 
lOtdyn at all scales r > 0.03r s . 

To simulate the effects of substructures on dark matter 
halos, we add ~ 1000 subhalos to the main halo. The 
subhalo masses are drawn from dn su b/dM su b cx M~V[, 
similar to those found in cosmological simulations [2j. 
Initially the subhalos are placed within the virial radius of 
the main halo either with a top-hat or r~ 2 radial number 
density distribution (see Fig. 1 for a comparison). The 
center-of-mass velocities of the subhalos are drawn from 
a local isotropic distribution identical to that of the CDM 
particles in the main halo. Simulations indicate that both 
CDM particles and subhalos are likely to develop mild 
velocity anisotropies in the outer part of the main halo 
[i"f)j | but this effect should be small. 

Point-Mass Subhalos. — In the simplest case, we repre- 
sent each subhalo as a point mass. This model is un- 
realistic and overestimates the relaxation effect since it 
ignores mass losses due to tidal stripping. However, it 
serves as a test case for the validity of the standard Chan- 
drasekhar formula (see Eq. [2] ) and approximates the ef- 
fects of dense baryonic clumps that can survive into halo 
centers. 

Fig. 1 shows the results of our point-mass subhalo sim- 
ulations in dimensionless units, which hold for halos of 
different masses at appropriately scaled cosmic times (see 
caption). It illustrates that point-mass subhalos can in- 
deed result in significant flattening in the inner p of a 
main halo within a few (inner) dynamical times. The 
amount of flattening is sensitive to the masses of the 
several most massive subhalos present in the main halo 
since these subhalos dominate the energy exchange with 
the dark matter in the main halo, as seen in Eq. (1). We 
have performed two identical runs - one having 999 point 
subhalos (dashed curves); the other having 996 point sub- 
halos (dotted) without the top 3 most massive subhalos 
in the 999 run - to test the effect of massive subhalos. 
The subhalos are placed initially within the main halo 
with r~ 2 distribution. Fig. 1 shows the heating of the 
main halo in the 996 run occurs at a later time (between 
2.8 and Q.9td yn ) and also leads to less flattening than 
the 999 run. The total subhalo mass in the 999 and 996 
runs is 7.02% and 3.29%M ma i„; the three most massive 
subhalos in the 999 run have masses 1.51%, 1.25%, and 
0.97%M mam . 

The two 999 subhalos runs in Fig. 1 illustrate the de- 
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FIG. 1: Evolution of the radial density profile of a main halo 
(with Cmain = 5.2) containing 996 vs. 999 point-mass subha- 
los. The dotted (996) and dashed (999) curves are for identical 
simulations except for the removal of the three most massive 
subhalos. The inner p(r) decreases with time from t = 
(solid) to 2.8t dyn (r a ) (upper of each pair) and 6.9td yn (r s ) 
(lower of each pair); td y n(r s ) is 0.072, 0.25, and 0.45 /i -1 Gyr 
for a 1O 8 M (at z = 4), 1O 12 M (z = 1), and 1O 14 M (z = 0) 
main halo. The dot-dashed curves compare the same 999 run 
as the dashed curves except the subhalo centers are placed 
initially with a tophat instead of r~ 2 distribution; the main 
halo here flattens later between 6.9 tdyn (upper dot-dashed) 
and 9.7 tdyn (lower). Without the subhalos, we have tested 
that the solid curve does not change over at least 10 tdyn- 



pendence of the relaxation timescales on the subhalo spa- 
tial distribution. In accordance with Eq. (1), the inner 
part of the main halo flattens more quickly for the r~ 2 
case than the top-hat case. For the latter, the initial 
main halo p(r) is unchanged through 6.9 td yn and then 
flattens quickly in three tdyn- Tne stabilized main halo 
profile (at 9.7 td yn \ bottom dot-dashed) is similar to the 
other 999 case, so the difference due to different subhalo 
spatial distributions is mainly in the timescales. 

Puffy Subhalos. — To model the subhalos more realisti- 
cally, we perform a series of simulations in which the 10 
most massive subhalos are given NFW profiles, while 989 
lower mass subhalos (all < 0.1% M ma i n ) are represented 
by point particles since they would suffer little mass loss. 
The total mass in subhalos in these cases is equal to ei- 
ther 7.02 or 10.3%M ma i„, where the ten most massive 
subhalos comprise 5.2% or 9.0%. The 7.02% model has 
the same subhalo mass spectrum as in the point-mass 
runs above, while the 10.3% model uses a different re- 
alization in which the three most massive subhalos are 
4.66, 2.09, and 1.00%M mam . 

Fig. 2 shows p(r) from our simulations of NFW subha- 
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FIG. 2: Radial density profile of a main halo (with Cmain = 
5.2) containing puffy subhalos. The panels compare p(r) of 
the main halo (left) vs main-plus-subhalos (right) initially 
(solid) and after 5.55 tdyn of evolution (other 3 curves). Three 
simulations with different subhalo concentration (c s ub = 15.6 
vs 31.2) and total subhalo masses (10.3% vs 7% of main halo 
mass) are shown. 



los with three combinations of subhalo concentration c su b 
and subhalo mass fraction. Since the subhalos can now 
shed mass in complicated ways, we compute p(r) both 
from the main halo particles only (left) and from all par- 
ticles (right). Flattening in the main halo p(r) is seen for 
all three cases. The amount of flattening is more severe 
when a higher mass fraction of the system is in subhalos, 
and when the subhalos have a higher c sub because more 
centrally concentrated subhalos suffer less tidal mass loss 
as they sink towards the main halo center. 

The inner cusp of the total mass, however, can steepen, 
remain the same, or flatten, depending on the competi- 
tion between the addition from subhalo masses deposited 
in the central regions and the removal of main halo par- 
ticles due to gravitational heating. The three models in 
Fig. 2 (right panel) showcase the three outcomes. The 
subhalos in the 7.02% model are not massive enough to 
add much mass, so both the main and total p(r) are 
flattened (~ r -o.75^ m ^ 6td yn - In contrast, in the 
model with 10.3% subhalo mass and c su b = 31.2, the 
mass added by the most massive subhalos (the top two 
have 4.66% and 2.09% M main ) more than compensate 
for the flattening in the main halo, leading to a steeper 
than r _1 inner cusp. Fig. 3 illustrates this model in more 
detail with five time outputs: most of the evolution oc- 
curs within one dynamical time after 2.2t £ ; y „ when tie 
most massive subhalos make their way to the center. |ll| 
showed that accreting one massive concentrated subhalo 
of l0%M ma i n can also produce a cusp in an initially cored 
halo. 
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FIG. 3: Time evolution of p of the main halo (left) and main- 
plus-subhalo (right) for the c„ u b = 31.2, 10.3% subhalo run 
in Fig. 2. Most of the evolution occurs between 2.22 and 
3.33 t dy 

n ■ 

We emphasize that the inner p(r) of a halo depends 
sensitively on the location of the halo center used to com- 
pute p{r). Although the initial momenta of the subha- 
los are drawn from an isotropic distribution, fluctuations 
typically introduce a small center-of-mass (COM) mo- 
tion for the entire system: a COM velocity ~ 2% of the 
main halo circular velocity was not uncommon, result- 
ing in COM offsets of ~ r s in ~ 10 td yn - Neglecting 
this effect and naively using a halo's initial center as the 
center for subsequent outputs would lead to a flattened 
p(r). We use a more physical halo center (e.g. iteratively 
determined COM, most bound particle, or COM of the 
500 most bound particles; all three give nearly identical 
results), which eliminates this spurious flattening. 

Timescale. — How do the timescales seen in the simu- 
lations compare with the simple energy exchange time 
predicted by Eq. QJ? The latter predicts 

\M t v\ _ 1 yj 

Irelax - | D(A £)| ~ ^2 l n A p b M b 2xe~^ ' 1 ' 

where the second equality assumes M t -C M b - In our 
study, this gives the timescale for heating the dark mat- 
ter particles from a background of subhalos of mass M su b, 
density p su b in the main halo, and COM velocity disper- 
sion a su b. Re-expressing it in terms of the main halo's 
virial mass M ma i n and radius r v , and 1-d velocity disper- 
sion at r v , cr(r v ), we obtain 

0.12 10 pcritMmai 

H(z)lnA psubMgub \a(r v ) J 2xe~ x2 ' 

where the Hubble time at redshift z is H^ 1 (z) = 
9.78Gyr/i- 1 [fi m (l + zf + Oy]- 1 / 2 . Let dn sub /dM sub oc 
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M s °^ be the subhalo mass function (assuming a > 1), f3 



be the ratio of the total mass in subhalos to M„ 



and 7 



be the ratio of the most massive subhalo to M main . We 
find p C ritM main /p sub M sub = (3 - a)/(2 - a)/(200/?7); 
for the 999 point-mass model shown in Fig. 1, this ra- 
tio is about 19. The local dynamical time at the scale 
and virial radius of an NFW halo with c ma in = 5.2 is 
tdyn(r s ) = 0-l4:t dyn (r v ) = 0.046 H^ 1 {z) (assuming unity 
for factors involving velocities, which is likely an under- 
estimate), so we find t re i ax ~ 2.3 H~ 1 (z) ~ 50td yn (r s ) ~ 
7tdyn{r v ). This is at least 5 times longer than the flat- 
tening timescale seen in the 999 point-mass top-hat simu- 
lation in Fig. 1. Eqs. (l)-(3), however, are valid only for 
a stationary, infinite, homogeneous background with a 
global Maxwellian velocity distribution || . In our study, 
the background is an ensemble of dark matter subhalos, 
themselves moving in a deeper main halo potential and 
experiencing dynamical friction and tidal mass losses. 
While Eqs. (l)-(3) elucidate the energy exchange between 
subhalos and dark matter particles, it is not surprising 
that they do not predict the exact timescales seen in sim- 
ulations. 

We have performed a test run with 1000 point-mass 
subhalos of equal mass where the total subhalo mass is 
15% M ma i n . This subhalo mass spectrum is unrealistic, 
but this run provides an additional test case and a com- 
parison case for a recent cluster galaxy study [12j . Since 
each subhalo mass in our test run is only 0.015% M ma j n , 
~ 100 times smaller than the most massive subhalos in 
Fig. 1, t re iax in Eq. (3) increases by a factor of ~ 100, 
too long to result in change in the inner halo profile. We 
indeed did not see any flattening over ~ 9td yn (at r s ) in 
our simulation. 

Implications. — Our series of controlled numerical ex- 
periments indicates that collisionless dark matter parti- 
cles in the inner parts of galaxy and cluster halos can gain 
energy through gravitational scatterings off concentrated 
dark matter subhalos, altering the inner density cusp of 
the main halo within a few dynamical times. These sub- 
halos appear ubiquitous in high resolution cosmological 
simulations and provide the source of fluctuations for the 
diffusion described by Eq. QJ. We have studied the evo- 
lution of halos under the influence of only one generation 
of subhalos, while real halos grow continuously by ac- 
cretion and mergers. The effects we have seen, however, 
suggest that fluctuations due to subhalos in parent ha- 
los are important for understanding the time evolution 
of dark matter density profiles and the halo-to-halo scat- 
ter of the inner cusp seen in recent ultra-high resolution 
cosmological simulations |l3|. We have shown that this 
scatter may be explained by subhalo accretion histories: 
when we allow for a population of subhalos of varying 



concentration and mass, the total inner profile of dark 
matter can either steepen or flatten. 

Recent observations of dwarf galaxy rotation curves 
based on CO and Ha emission find significant variations 
in the inner profile, ranging from a core to ~ r _1 |l4j |. 
While baryonic physics can influence the central mass 
profile, the purely gravitational physics of subhalo scat- 
tering studied here may also accommodate the variations 
seen in these observations. Conversely, maintaining a 
stable and universal inner profile would require a "quies- 
cent" accretion history not involving concentrated mas- 
sive subhalos. Halos in cosmological models with trun- 
cated small-scale pow er ( e.g. warm dark matter) contain 
many fewer satellites ; their inner cusps should there- 
fore be less prone to variations due to subhalo accretion. 
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